Mind-Wandering

Model Fitting

Data are loaded from .Rdata.

setwd(getwd())
load("tms_data_preprocessed.Rdata")

bname = "tms_analyses"

models_task <- list(
  formula(probe.response ~ zbv * zlog.apen + (1|subj/condition)),
  formula(probe.response ~ zbv * zlog.apen + probeix + (1|subj/condition)),
  formula(probe.response ~ zbv * zlog.apen + probeix + block_num + (1|subj/condition)),
  formula(probe.response ~ zbv * zlog.apen + probeix + block_num + condition + (1|subj/condition)),
  formula(probe.response ~ zbv * zlog.apen + probeix + block_num + condition + randomization + (1|subj/condition))
)


descriptions_task=c(
  "BV x AE", 
  "BV x AE + trial", 
  "BV x AE + trial + block", 
  "BV x AE + trial + block + condition",
  "BV x AE + trial + block + condition + randomization"
  )

Model Diagnostics

loos_task <- load.cache.var("loos_task", base = bname) # load loos
CACHE> loading loos_task from cache/vars/tms_analyses_loos_task.RData
mod_selection_res <- as.data.frame(loos_task$diffs) %>% select(c(elpd_diff, se_diff)) # from loos select elpd_diff and se_diff, convert to data frame
models_elpd_diff <- c(
  "BV x AE + probe + block",
  "BV x AE + probe + block + condition",
  "BV x AE + probe + block + condition + randomization",
  "BV x AE + probe",
  "BV x AE"
  ) #List models in the same order as in elpd_diff. Those are labels for the graphical table
mod_selection_res <- data.frame(models_elpd_diff, mod_selection_res$elpd_diff, mod_selection_res$se_diff) # create a data frame with model labels and loo diagnostics
names(mod_selection_res) <- c("Model", "elpd_diff", "se_diff") # name columns appropriately

mod_selection_res = mod_selection_res  %>% 
  mutate(across(where(is.numeric), ~comma(., accuracy=0.01))) %>% 
  rename_all(~gsub("\\.", " ", .)) # format table

alignment= map_chr(mod_selection_res, ~ifelse(class(.x)=="numeric", "r","l")) # align rows
mod_selection_res %>%  # plot model selection table
  kbl(caption="LOO-CV: MW", align=alignment) %>% 
  kable_classic(full_width=TRUE, html_font="Times") %>% 
  row_spec(2, bold = TRUE, color = "red")
LOO-CV: MW
Model elpd_diff se_diff
BV x AE + probe + block 0.00 0.00
BV x AE + probe + block + condition -0.24 1.25
BV x AE + probe + block + condition + randomization -0.29 1.26
BV x AE + probe -13.58 5.44
BV x AE -17.27 6.33

Model Coefficients: lots of uncertainty around stimulation (95%-intervals)

color_scheme_set("viridisE")
bayesplot_theme_set(theme_bw(base_size = 15, base_family = "sans"))

mod_task03 <- load.cache.var("mod_task03", base = bname)

CACHE> loading mod_task03 from cache/vars/tms_analyses_mod_task03.RData

mcmc_intervals(mod_task03, prob_outer = 0.95, pars = c("b_zbv",
                         "b_zlog.apen",
                         "b_block_num",
                         "b_probeix",
                         "b_zbv:zlog.apen",
                        # "b_visit2",
                         "b_conditionsham_rhTMS",
                         "b_conditionsham_arrhTMS",
                         "b_conditionactive_rhTMS",
                         "b_conditionactive_arrhTMS"
                         )) +
  ggplot2::scale_y_discrete(labels = c("b_zbv" = "BV",
                                       "b_zlog.apen" = "AE",
                                       "b_block_num" = "Block",
                                       "b_probeix" = "Probe number",
                                       "b_zbv:zlog.apen" = "BV x AE",
                                      # "b_visit2" = "Visit",
                                       "b_conditionsham_rhTMS" = "Sham rhTMS",
                                       "b_conditionsham_arrhTMS" = "Sham arrhTMS",
                                       "b_conditionactive_rhTMS" = "Active rhTMS",
                                       "b_conditionactive_arrhTMS" = "Active arrhTMS")) + 
geom_vline(xintercept = 0,color = "red", size=0.75)  + theme(axis.text.y = element_text(size=18, face = "plain", colour = c('black','black', 'black', 'black', 'black', 'black', 'black', 'red', 'black'))) + ggtitle("Model Coefficients: MW")
Scale for 'y' is already present. Adding another scale for 'y', which will
replace the existing scale.
Fig. 1. Model 7: On-task Score: coefficients for all parameters of interest.

Fig. 1. Model 7: On-task Score: coefficients for all parameters of interest. 

mw_hyp = hypothesis(mod_task03, c("zbv<0", "zlog.apen>0", "block_num<0", "probeix<0", "zbv:zlog.apen<0", "conditionactive_rhTMS>0", "conditionsham_rhTMS>0", "conditionactive_arrhTMS>0", "conditionsham_arrhTMS>0"), alpha = 0.025) # test hypotheses for model 03
mw_hyp_table <- mw_hyp$hypothesis %>% select(-Star) # select all columns from the resulting object except for Star
  
alignment= map_chr(mw_hyp_table, ~ifelse(class(.x)=="numeric", "r","l")) # alogn rows
hypotheses <- c("BV < 0", "AE > 0", "Block < 0", "Probe number < 0", "BV x AE< 0", "Active rhTMS > 0", "Sham rhTMS > 0", "Active arrhTMS > 0", "Sham arrhTMS > 0") # labels for the hypothesis column
mw_hyp_table = mw_hyp_table  %>% # table formatting
  mutate(across(where(is.numeric), ~comma(., accuracy=0.01)), Hypothesis = hypotheses) %>% 
  rename_all(~gsub("\\.", " ", .)) 
mw_hyp_table %>% # plot table
  kbl(caption="Results: MW", align=alignment) %>% 
  kable_classic(full_width=TRUE, html_font="Times") %>% 
  row_spec(c(7, 8, 9), bold = TRUE) %>% 
  row_spec(6, bold = TRUE, color = "red")
Results: MW
Hypothesis Estimate Est Error CI Lower CI Upper Evid Ratio Post Prob
BV < 0 -0.03 0.05 -0.13 0.06 3.37 0.77
AE > 0 0.06 0.05 -0.03 0.15 9.58 0.91
Block < 0 -0.11 0.02 -0.15 -0.07 Inf 1.00
Probe number < 0 -0.04 0.01 -0.07 -0.02 713.29 1.00
BV x AE< 0 -0.04 0.05 -0.14 0.06 3.75 0.79
Active rhTMS > 0 0.28 0.26 -0.23 0.79 6.52 0.87
Sham rhTMS > 0 -0.23 0.26 -0.75 0.26 0.22 0.18
Active arrhTMS > 0 0.06 0.26 -0.44 0.58 1.43 0.59
Sham arrhTMS > 0 -0.32 0.26 -0.84 0.19 0.12 0.11

Long data frame for plotting

tms_data.nback_z <- tms_data.nback %>%
  reshape2::melt(id.vars = c( "focus", "probe.response", "intention", "somnolence", "condition", "visit", "subj"), measure.vars = c("zlog.apen", "zbv"), varnames = c("Variable", "Score"))

data.table::setnames(tms_data.nback_z, old = c("focus", "probe.response", "intention", "somnolence", "condition", "visit", "subj",'variable', "value"), new = c('Focus', "On-task Score", "Intention", "Alertness", "Condition", "Visit", "Subject",'Measure', "Z-score"))

On-task Score Across Conditons

Fig. 2. MW Mean ± 1SE across conditions.

Fig. 2. MW Mean ± 1SE across conditions.

Model fitting: BV

Model Diagnostics

loos_bv=load.cache.var("loos_bv", base=bname)
CACHE> loading loos_bv from cache/vars/tms_analyses_loos_bv.RData
# same procedure for BV and AE
as.data.frame(loos_bv$diffs)
mod_selection_res <- as.data.frame(loos_bv$diffs) %>% select(c(elpd_diff, se_diff))
models_elpd_diff <- c(
  "probe + block + condition + visit",
  "probe + block",
  "probe + block + condition",
  "probe"
  )
mod_selection_res <- data.frame(models_elpd_diff, mod_selection_res$elpd_diff, mod_selection_res$se_diff)
names(mod_selection_res) <- c("Model", "elpd_diff", "se_diff")

mod_selection_res = mod_selection_res  %>% 
  mutate(across(where(is.numeric), ~comma(., accuracy=0.01))) %>% 
  rename_all(~gsub("\\.", " ", .)) 

alignment= map_chr(mod_selection_res, ~ifelse(class(.x)=="numeric", "r","l"))
         elpd_diff  se_diff  elpd_loo se_elpd_loo    p_loo se_p_loo    looic
mod_bv03  0.000000 0.000000 -749.5354    37.47914 44.65426 1.245530 1499.071
mod_bv01 -2.226843 4.354602 -751.7623    37.03687 40.26908 1.058170 1503.525
mod_bv02 -3.282269 3.992372 -752.8177    36.94252 41.68978 1.117909 1505.635
mod_bv00 -4.509851 5.198430 -754.0453    37.10229 38.84448 1.067001 1508.091
         se_looic
mod_bv03 74.95829
mod_bv01 74.07375
mod_bv02 73.88504
mod_bv00 74.20458
mod_selection_res %>% 
  kbl(caption="LOO-CV: BV", align=alignment) %>% 
  kable_classic(full_width=TRUE, html_font="Times") %>% 
  row_spec(1, bold = TRUE, color = "red")
LOO-CV: BV
Model elpd_diff se_diff
probe + block + condition + visit 0.00 0.00
probe + block -2.23 4.35
probe + block + condition -3.28 3.99
probe -4.51 5.20

BV: Model Coefficients

color_scheme_set("viridisE")
bayesplot_theme_set(theme_bw(base_size = 15, base_family = "sans"))

mod_bv03 <- load.cache.var("mod_bv03", base = bname)
mcmc_intervals(mod_bv03, prob_outer = 0.95, pars = c(
                         "b_block_num",
                         "b_probeix",
                          "b_visit2",
                         "b_conditionsham_rhTMS",
                         "b_conditionsham_arrhTMS",
                         "b_conditionactive_rhTMS",
                         "b_conditionactive_arrhTMS"
                         )) +
  ggplot2::scale_y_discrete(labels = c("b_block_num" = "Block",
                                       "b_probeix" = "Probe number",
                                       "b_visit2" = "Visit",
                                       "b_conditionsham_rhTMS" = "Sham rhTMS",
                                       "b_conditionsham_arrhTMS" = "Sham arrhTMS",
                                       "b_conditionactive_rhTMS" = "Active rhTMS",
                                       "b_conditionactive_arrhTMS" = "Active arrhTMS")) + 
geom_vline(xintercept = 0,color = "red", size=0.75)  + theme(axis.text.y = element_text(size=18, face = "plain", colour = c('black', 'black', 'black', 'black', 'black', 'red', 'black'))) + ggtitle("Model Coefficients: BV")
Scale for 'y' is already present. Adding another scale for 'y', which will
replace the existing scale.
Fig. 6. BV: coefficients for all parameters of interest.

Fig. 6. BV: coefficients for all parameters of interest. 

CACHE> loading mod_bv03 from cache/vars/tms_analyses_mod_bv03.RData
bv_hyp = hypothesis(mod_bv03, c("block_num>0", "probeix>0", "conditionactive_rhTMS>0", "conditionsham_rhTMS>0", "conditionactive_arrhTMS>0", "conditionsham_arrhTMS>0"), alpha = 0.025)
bv_hyp_table <- bv_hyp$hypothesis %>% select(-Star)
  
alignment= map_chr(bv_hyp_table, ~ifelse(class(.x)=="numeric", "r","l"))
hypotheses <- c("Block > 0", "Probe number > 0", "Active rhTMS > 0", "Sham rhTMS > 0", "Active arrhTMS > 0", "Sham arrhTMS > 0")
bv_hyp_table = bv_hyp_table  %>% 
  mutate(across(where(is.numeric), ~comma(., accuracy=0.01)), Hypothesis = hypotheses) %>% 
  rename_all(~gsub("\\.", " ", .)) 
bv_hyp_table %>% 
  kbl(caption="Results: BV", align=alignment) %>% 
  kable_classic(full_width=TRUE, html_font="Times") %>% 
  row_spec(c(4, 5, 6), bold = TRUE) %>% 
  row_spec(3, bold = TRUE, color = "red")
Results: BV
Hypothesis Estimate Est Error CI Lower CI Upper Evid Ratio Post Prob
Block > 0 -0.02 0.01 -0.04 0.00 0.01 0.01
Probe number > 0 0.00 0.01 -0.01 0.01 1.17 0.54
Active rhTMS > 0 0.18 0.09 0.01 0.35 50.02 0.98
Sham rhTMS > 0 0.04 0.09 -0.12 0.22 2.36 0.70
Active arrhTMS > 0 0.13 0.09 -0.04 0.30 13.90 0.93
Sham arrhTMS > 0 0.04 0.08 -0.12 0.21 2.21 0.69

Model fitting: AE

Model Diagnostics

loos_apen=load.cache.var("loos_apen", base=bname)
CACHE> loading loos_apen from cache/vars/tms_analyses_loos_apen.RData

Model Coefficients: lots of uncertainty around stimulation (90%-intervals)

mod_selection_res <- as.data.frame(loos_apen$diffs) %>% select(c(elpd_diff, se_diff))
models_elpd_diff <- c(
  "probe + block",
  "probe",
  "probe + block + condition",
  "probe + block + condition + visit"
  )
mod_selection_res <- data.frame(models_elpd_diff, mod_selection_res$elpd_diff, mod_selection_res$se_diff)
names(mod_selection_res) <- c("Model", "elpd_diff", "se_diff")

mod_selection_res = mod_selection_res  %>% 
  mutate(across(where(is.numeric), ~comma(., accuracy=0.01))) %>% 
  rename_all(~gsub("\\.", " ", .)) 

alignment= map_chr(mod_selection_res, ~ifelse(class(.x)=="numeric", "r","l"))
mod_selection_res %>% 
  kbl(caption="LOO-CV: AE", align=alignment) %>% 
  kable_classic(full_width=TRUE, html_font="Times") %>% 
  row_spec(4, bold = TRUE, color = "red")
LOO-CV: AE
Model elpd_diff se_diff
probe + block 0.00 0.00
probe -0.77 1.99
probe + block + condition -0.97 1.97
probe + block + condition + visit -1.45 2.01 
color_scheme_set("viridisE")
bayesplot_theme_set(theme_bw(base_size = 15, base_family = "sans"))

mod_apen03 <- load.cache.var("mod_apen03", base = bname)
mcmc_intervals(mod_apen03, prob_outer = 0.95, pars = c(
                         "b_block_num",
                         "b_probeix",
                          "b_visit2",
                         "b_conditionsham_rhTMS",
                         "b_conditionsham_arrhTMS",
                         "b_conditionactive_rhTMS",
                         "b_conditionactive_arrhTMS"
                         )) +
  ggplot2::scale_y_discrete(labels = c("b_block_num" = "Block",
                                       "b_probeix" = "Probe number",
                                       "b_visit2" = "Visit",
                                       "b_conditionsham_rhTMS" = "Sham rhTMS",
                                       "b_conditionsham_arrhTMS" = "Sham arrhTMS",
                                       "b_conditionactive_rhTMS" = "Active rhTMS",
                                       "b_conditionactive_arrhTMS" = "Active arrhTMS")) + 
geom_vline(xintercept = 0,color = "red", size=0.75)  + theme(axis.text.y = element_text(size=18, face = "plain", colour = c('black', 'black', 'black', 'black', 'black', 'red', 'black'))) + ggtitle("Model Coefficients: AE")
Scale for 'y' is already present. Adding another scale for 'y', which will
replace the existing scale.
Warning: Vectorized input to `element_text()` is not officially supported.
Results may be unexpected or may change in future versions of ggplot2.
CACHE> loading mod_apen03 from cache/vars/tms_analyses_mod_apen03.RData
Fig. 7. AE: coefficients for all parameters of interest.

Fig. 7. AE: coefficients for all parameters of interest. 

apen_hyp = hypothesis(mod_apen03, c("block_num<0", "probeix<0", "conditionactive_rhTMS<0", "conditionsham_rhTMS<0", "conditionactive_arrhTMS<0", "conditionsham_arrhTMS<0"), alpha = 0.025)
apen_hyp_table <- apen_hyp$hypothesis %>% select(-Star)
  
alignment= map_chr(apen_hyp_table, ~ifelse(class(.x)=="numeric", "r","l"))
hypotheses <- c("Block < 0", "Probe number < 0", "Active rhTMS < 0", "Sham rhTMS < 0", "Active arrhTMS < 0", "Sham arrhTMS < 0")
apen_hyp_table = apen_hyp_table  %>% 
  mutate(across(where(is.numeric), ~comma(., accuracy=0.01)), Hypothesis = hypotheses) %>% 
  rename_all(~gsub("\\.", " ", .)) 
apen_hyp_table %>% 
  kbl(caption="Results: AE", align=alignment) %>% 
  kable_classic(full_width=TRUE, html_font="Times") %>% 
  row_spec(c(4, 5, 6), bold = TRUE) %>% 
  row_spec(3, bold = TRUE, color = "red")
Results: AE
Hypothesis Estimate Est Error CI Lower CI Upper Evid Ratio Post Prob
Block < 0 -0.02 0.01 -0.05 0.00 23.81 0.96
Probe number < 0 0.02 0.01 0.00 0.04 0.01 0.01
Active rhTMS < 0 -0.27 0.12 -0.50 -0.04 83.03 0.99
Sham rhTMS < 0 -0.22 0.12 -0.45 0.01 32.78 0.97
Active arrhTMS < 0 -0.23 0.12 -0.47 -0.01 43.25 0.98
Sham arrhTMS < 0 -0.12 0.12 -0.35 0.11 5.97 0.86

Non-Parametric ANOVA

Shapiro Test

# Check data normality assumption. Spoiler alert: violated.

shapiro.test(tms_data.nback$probe.response)
shapiro.test(tms_data.nback$zlog.apen)
shapiro.test(tms_data.nback$bv)

    Shapiro-Wilk normality test

data:  tms_data.nback$probe.response
W = 0.87298, p-value < 2.2e-16


    Shapiro-Wilk normality test

data:  tms_data.nback$zlog.apen
W = 0.9406, p-value < 2.2e-16


    Shapiro-Wilk normality test

data:  tms_data.nback$bv
W = 0.5741, p-value < 2.2e-16

Distributions are non-normal.

MW: Kruskal-Wallis Test & Post-Hoc Multiple Comparisons

kruskal.test(probe.response ~ condition, data = tms_data.nback) # groups are different

    Kruskal-Wallis rank sum test

data:  probe.response by condition
Kruskal-Wallis chi-squared = 15.76, df = 4, p-value = 0.003359

MW: Wilcoxon tests to zoom in

tms_data.nback$condition <- fct_relevel(tms_data.nback$condition, "baseline", "active_rhTMS", "active_arrhTMS", "sham_rhTMS", "sham_arrhTMS") # relevel condition factors

tms_comparison_list <- list(c("baseline", "active_rhTMS"), c("baseline", "active_arrhTMS"), c("baseline", "sham_rhTMS"), c("baseline", "sham_arrhTMS"), c("active_rhTMS", "active_arrhTMS"), c("active_rhTMS", "sham_rhTMS"), c("active_rhTMS", "sham_arrhTMS"), c("active_arrhTMS", "sham_rhTMS"), c("active_arrhTMS", "sham_arrhTMS"), c("sham_rhTMS", "sham_arrhTMS")) # define the exact comparison pairs

pairwise_wilcox_results <- data.frame(matrix(ncol = 4, nrow = 10)) # create an empty data frame for Wilcoxon results
colnames(pairwise_wilcox_results) <- c("Comparison", "W-statistic", "p-value", "Effect size, r") # name the columns
pairwise_wilcox_results$Comparison <- tms_comparison_list # assign comparisons to the corresponding column
eff_sizes_mw <- rstatix::wilcox_effsize(tms_data.nback, probe.response ~ condition, paired = FALSE) # compute effect size estimates

i = 1
 for (comparison in tms_comparison_list) { # perform Wilcox test for each comparison pair
   r <- wilcox.test(filter(tms_data.nback, condition == comparison[1])$probe.response, filter(tms_data.nback, condition == comparison[2])$probe.response, paired = FALSE)
  pairwise_wilcox_results[i, 2:4] <- rbind(bind_cols(r$statistic, r$p.value, eff_sizes_mw$effsize[[i]])) #bind rows
  i = i + 1
  
 }

comparisons <- c("baseline vs. active rhTMS", "baseline vs. active arrhTMS", "baseline vs. sham rhTMS", "baseline vs. sham arrhTMS", "active rhTMS vs. active arrhTMS", "active rhTMS vs. sham rhTMS", "active rhTMS vs. sham arrhTMS", "active arrhTMS vs. sham rhTMS", "active arrhTMS vs. sham arrhTMS", "sham rhTMS vs. sham arrhTMS") # labels for comparisons

pairwise_wilcox_results <- pairwise_wilcox_results %>% # format table
  mutate(across("p-value", ~comma(., accuracy=0.01)), across("Effect size, r", ~comma(., accuracy=0.01)),Comparison = comparisons) %>% 
  rename_all(~gsub("\\.", " ", .)) 
pairwise_wilcox_results %>% # plot table
  kbl(caption="Pairwise Wilcoxon test: MW", align=alignment) %>% 
  kable_classic(full_width=TRUE, html_font="Times") %>% 
  row_spec(c(3, 4, 6, 7), bold = TRUE) %>% 
  row_spec(1, color = "red")
Pairwise Wilcoxon test: MW
Comparison W-statistic p-value Effect size, r
baseline vs. active rhTMS 18312 0.41 0.04
baseline vs. active arrhTMS 20483 0.23 0.06
baseline vs. sham rhTMS 21615 0.02 0.11
baseline vs. sham arrhTMS 22196 0.00 0.14
active rhTMS vs. active arrhTMS 14191 0.07 0.10
active rhTMS vs. sham rhTMS 14899 0.01 0.15
active rhTMS vs. sham arrhTMS 15308 0.00 0.18
active arrhTMS vs. sham rhTMS 13626 0.29 0.06
active arrhTMS vs. sham arrhTMS 13944 0.14 0.08
sham rhTMS vs. sham arrhTMS 13088 0.71 0.02

AE: Kruskal-Wallis Test & Post-Hoc Multiple Comparisons

# do the same for BV and AE
kruskal.test(zlog.apen ~ condition, data = tms_data.nback) # groups are different

    Kruskal-Wallis rank sum test

data:  zlog.apen by condition
Kruskal-Wallis chi-squared = 9.6253, df = 4, p-value = 0.04724
pairwise_wilcox_results <- data.frame(matrix(ncol = 4, nrow = 10))
colnames(pairwise_wilcox_results) <- c("Comparison", "W-statistic", "p-value", "Effect size, r")
pairwise_wilcox_results$Comparison <- tms_comparison_list
eff_sizes_mw <- rstatix::wilcox_effsize(tms_data.nback, zlog.apen ~ condition, paired = FALSE)
i = 1
 for (comparison in tms_comparison_list) {
   r <- wilcox.test(filter(tms_data.nback, condition == comparison[1])$zlog.apen, filter(tms_data.nback, condition == comparison[2])$zlog.apen, paired = FALSE)
  pairwise_wilcox_results[i, 2:4] <- rbind(bind_cols(r$statistic, r$p.value, eff_sizes_mw$effsize[[i]]))
  i = i + 1
  
 }

pairwise_wilcox_results <- pairwise_wilcox_results %>% 
  mutate(across("p-value", ~comma(., accuracy=0.01)), across("Effect size, r", ~comma(., accuracy=0.01)),Comparison = comparisons) %>% 
  rename_all(~gsub("\\.", " ", .)) 
pairwise_wilcox_results %>% 
  kbl(caption="Pairwise Wilcoxon test: AE", align=alignment) %>% 
  kable_classic(full_width=TRUE, html_font="Times") %>% 
  row_spec(c(1, 2), bold = TRUE) %>% 
  row_spec(1, color = "red")
Pairwise Wilcoxon test: AE
Comparison W-statistic p-value Effect size, r
baseline vs. active rhTMS 22138.5 0.01 0.13
baseline vs. active arrhTMS 21689.0 0.03 0.11
baseline vs. sham rhTMS 21311.5 0.06 0.09
baseline vs. sham arrhTMS 19889.5 0.54 0.03
active rhTMS vs. active arrhTMS 12561.5 0.77 0.02
active rhTMS vs. sham rhTMS 12431.0 0.66 0.02
active rhTMS vs. sham arrhTMS 11337.5 0.08 0.10
active arrhTMS vs. sham rhTMS 12673.5 0.88 0.01
active arrhTMS vs. sham arrhTMS 11567.5 0.14 0.08
sham rhTMS vs. sham arrhTMS 11767.5 0.21 0.07

BV: Kruskal-Wallis Test & Post-Hoc Multiple Comparisons

pairwise_wilcox_results <- data.frame(matrix(ncol = 3, nrow = 10))
colnames(pairwise_wilcox_results) <- c("Comparison", "W-statistic", "p-value")
pairwise_wilcox_results$Comparison <- tms_comparison_list

i = 1
 for (comparison in tms_comparison_list) {
   r <- wilcox.test(filter(tms_data.nback, condition == comparison[1])$zbv, filter(tms_data.nback, condition == comparison[2])$zbv, paired = FALSE)
  pairwise_wilcox_results[i, 2:3] <- rbind( bind_cols(r$statistic, r$p.value))
  i = i + 1
  
 }

pairwise_wilcox_results <- pairwise_wilcox_results %>% 
  mutate(across("p-value", ~comma(., accuracy=0.01)), Comparison = comparisons) %>% 
  rename_all(~gsub("\\.", " ", .)) 
pairwise_wilcox_results %>% 
  kbl(caption="Pairwise Wilcoxon test: BV", align=alignment) %>% 
  kable_classic(full_width=TRUE, html_font="Times") %>% 
  row_spec(c(2), bold = TRUE) %>% 
  row_spec(1, color = "red")
Pairwise Wilcoxon test: BV
Comparison W-statistic p-value
baseline vs. active rhTMS 16467 0.02
baseline vs. active arrhTMS 16813 0.04
baseline vs. sham rhTMS 18209 0.38
baseline vs. sham arrhTMS 17964 0.28
active rhTMS vs. active arrhTMS 13154 0.67
active rhTMS vs. sham rhTMS 13859 0.20
active rhTMS vs. sham arrhTMS 13786 0.23
active arrhTMS vs. sham rhTMS 13625 0.32
active arrhTMS vs. sham arrhTMS 13615 0.32
sham rhTMS vs. sham arrhTMS 12678 0.88

AE: Wilcoxon test to zoom in at the significant contrast

wilcox.test(filter(tms_data.nback_z, (`Condition` == "baseline" & `Measure` == "zlog.apen"))$`Z-score`,
            filter(tms_data.nback_z, (`Condition` == "active_rhTMS" & `Measure` == "zlog.apen"))$`Z-score`,
            paired = FALSE)

rstatix::wilcox_effsize(tms_data.nback, zlog.apen ~ condition, paired = FALSE)

    Wilcoxon rank sum test with continuity correction

data:  filter(tms_data.nback_z, (Condition == "baseline" & Measure == "zlog.apen"))$`Z-score` and filter(tms_data.nback_z, (Condition == "active_rhTMS" & Measure == "zlog.apen"))$`Z-score`
W = 22138, p-value = 0.009497
alternative hypothesis: true location shift is not equal to 0

# A tibble: 10 x 7
   .y.       group1         group2         effsize    n1    n2 magnitude
 * <chr>     <chr>          <chr>            <dbl> <int> <int> <ord>    
 1 zlog.apen baseline       active_rhTMS   0.130     240   160 small    
 2 zlog.apen baseline       active_arrhTMS 0.110     240   160 small    
 3 zlog.apen baseline       sham_rhTMS     0.0932    240   160 small    
 4 zlog.apen baseline       sham_arrhTMS   0.0304    240   160 small    
 5 zlog.apen active_rhTMS   active_arrhTMS 0.0161    160   160 small    
 6 zlog.apen active_rhTMS   sham_rhTMS     0.0249    160   160 small    
 7 zlog.apen active_rhTMS   sham_arrhTMS   0.0988    160   160 small    
 8 zlog.apen active_arrhTMS sham_rhTMS     0.00855   160   160 small    
 9 zlog.apen active_arrhTMS sham_arrhTMS   0.0833    160   160 small    
10 zlog.apen sham_rhTMS     sham_arrhTMS   0.0697    160   160 small

BV: Kruskal-Wallis Test & Post-Hoc Multiple Comparisons

kruskal.test(zbv ~ condition, data = tms_data.nback) # groups are NOT different

    Kruskal-Wallis rank sum test

data:  zbv by condition
Kruskal-Wallis chi-squared = 7.4528, df = 4, p-value = 0.1138

A Closer Look at Task Performance

tms_data.nback_z %>%
  ggplot(aes(x= `Condition`, y =`Z-score`, group = `Measure`, shape = `Condition`, color = `Measure`)) + 
  geom_pointrange(stat="summary", fun.data=mean_se, fun.args = list(mult=1), position=position_dodge(0.05)) +
  geom_line(stat="summary", fun.data=mean_se, fun.args = list(mult=2)) +
  scale_color_manual(labels = c("AE", "BV"), values = c("blue", "red")) +
  theme_set(theme_bw()) +
  theme(text = element_text(size = 15)) +
  theme(axis.text.x = element_text(angle = 20,  vjust = 0.5, hjust=0.4))
Fig. 3. Mean ± 1SE changes of AE and BV across conditions.

Fig. 3. Mean ± 1SE changes of AE and BV across conditions. 

tms_data.nback_z %>%
  ggplot(aes(x= `On-task Score`, y =`Z-score`, group = `Measure`, color = `Measure`)) + 
  geom_pointrange(stat="summary", fun.data=mean_se, fun.args = list(mult=1), position=position_dodge(0.05)) +
  geom_line(stat="summary", fun.data=mean_se, fun.args = list(mult=2)) +
  facet_wrap(~`Condition`) +
  scale_color_manual(labels = c("AE", "BV"), values = c("blue", "red")) +
  theme(text = element_text(size = 12))
Fig. 4. AE & BV versus On-task score: Mean ± 1SE across conditions. Our subjects are either bad at tapping with the metronome or misunderstood this part of the task.

Fig. 4. AE & BV versus On-task score: Mean ± 1SE across conditions. Our subjects are either bad at tapping with the metronome or misunderstood this part of the task.

AE x BV Interaction Across Conditions

tms_data.nback_z %>%
  ggplot(aes(x= `Focus`, y =`Z-score`,  group = `Measure`, color=`Measure`)) + 
  geom_pointrange(aes(shape = `Condition`),stat="summary", fun.data=mean_se, fun.args = list(mult=1), position=position_dodge(0.05)) +
  geom_line(stat="summary", fun.data=mean_se, fun.args = list(mult=2)) +
  scale_color_manual(labels = c("AE", "BV"), values = c("blue", "red")) +
  facet_wrap(~ `Condition`)
Fig. 5. AE & BV versus On-task score: BV x AE interaction versus on-off task states. acrive rhTMS seems to flip the interaction: subjects have worse performance when they think they are on task and vice-versa.

Fig. 5. AE & BV versus On-task score: BV x AE interaction versus on-off task states. acrive rhTMS seems to flip the interaction: subjects have worse performance when they think they are on task and vice-versa.

Appendices

Histograms

tms_data.nback %>%
  ggplot(aes(x = zbv)) +
  geom_histogram(binwidth = 0.1) +
  labs(y = "Count", x = "Z-score: BV") +
  theme(text = element_text(size = 18)) 

tms_data.nback %>%
  ggplot(aes(x = zlog.apen)) +
  geom_histogram(binwidth = 0.1) +
  labs(y = "Count", x = "Z-score: AE")  +
  theme(text = element_text(size = 18)) 

tms_data.nback %>%
  ggplot(aes(x = probe.response)) +
  geom_histogram(binwidth = 0.5) +
  labs(y = "Count", x = "On-task score")  +
  theme(text = element_text(size = 18))
AE and BV distributionsAE and BV distributionsAE and BV distributions

AE and BV distributions

Subject-wise. Better zoom in RStudio output

tms_data.nback_z %>%
  ggplot(aes(x= `Condition`, y =`On-task Score`, group = 1)) + 
  geom_pointrange(stat="summary", fun.data=mean_se, fun.args = list(mult=1), position=position_dodge(0.05)) +
  geom_line(stat="summary", fun.data=mean_se, fun.args = list(mult=2)) +
  facet_wrap(~ `Subject`, scales = "free", ncol = 2, as.table = FALSE) +
  theme_set(theme_bw()) +
  geom_vline(xintercept = "active_rhTMS", color = "red", size=0.5, linetype= "dotdash") +
  theme(axis.text.x = element_text(angle = 20,  vjust = 0.5, hjust=0.4)) +
  theme(plot.margin = margin(1,1,1.5,1.2, "cm"))
AE and BV distributions

AE and BV distributions 

tms_data.nback_z %>%
  ggplot(aes(x= `On-task Score`, y =`Z-score`, group = `Measure`, color = `Measure`)) + 
  geom_pointrange(stat="summary", fun.data=mean_se, fun.args = list(mult=1), position=position_dodge(0.05)) +
  geom_line(stat="summary", fun.data=mean_se, fun.args = list(mult=2)) +
  facet_wrap(~`Condition`) +
  scale_color_manual(labels = c("AE", "BV"), values = c("blue", "red")) +
  theme(text = element_text(size = 15)) 

tms_data.nback_z %>%
  ggplot(aes(x= `Focus`, y =`Z-score`,  group = `Measure`, color=`Measure`)) + 
  geom_pointrange(stat="summary", fun.data=mean_se, fun.args = list(mult=1), position=position_dodge(0.05)) +
  geom_line(stat="summary", fun.data=mean_se, fun.args = list(mult=2)) +
  scale_color_manual(labels = c("AE", "BV"), values = c("blue", "red")) +
  facet_wrap(~ `Condition`)+
  theme(text = element_text(size = 18))